# Cite Script from:
# De Moraes, R. M., Istoe, R. S. C., Miranda, V. A. (2023). How political skill and work engagement differ by hierarchical level: evidence from Brazil. Revista Brasileira de Gestão de Negócios, 25(3), firstpage-lastpage.

# Packages:

if(!require(dplyr)) install.packages("dplyr")
library(dplyr)                                
if(!require(car)) install.packages("car")   
library(car)                                
if(!require(rstatix)) install.packages("rstatix") 
library(rstatix)                                
if(!require(emmeans)) install.packages("emmeans") 
library(emmeans)
if(!require(ggplot2)) install.packages("ggplot2") 
library(ggplot2)
if(!require(MVN)) install.packages("MVN") 
library(MVN)
if(!require(GGally)) install.packages("GGally") 
library(GGally)
if (!require(boot)) install.packages("boot")
library(boot)

# load Dataset

dados <- read.csv2('DadosR.csv', stringsAsFactors = T,
                   fileEncoding = "latin1")  
View(dados)                                  
glimpse(dados)                               
head(dados, 10)

# Sample Charac.

media_age <- mean(dados$Age)
desvio_padrao_age <- sd(dados$Age)

media_jobten <- mean(dados$JobTen)
desvio_padrao_jobten <- sd(dados$JobTen)

cat("1) Média e desvio padrão para Age:\n")
cat("Média:", media_age, "\n")
cat("Desvio Padrão:", desvio_padrao_age, "\n\n")

cat("Média e desvio padrão para JobTen:\n")
cat("Média:", media_jobten, "\n")
cat("Desvio Padrão:", desvio_padrao_jobten, "\n\n")

# 2) Freq. and % Gender, Idade_G, Educ, Atividades e NivH
frequencia_gender <- table(dados$Gender)
percentual_gender <- prop.table(frequencia_gender) * 100

frequencia_idade_g <- table(dados$Idade_G)
percentual_idade_g <- prop.table(frequencia_idade_g) * 100

frequencia_educ <- table(dados$Educ)
percentual_educ <- prop.table(frequencia_educ) * 100

frequencia_atividades <- table(dados$Atividades)
percentual_atividades <- prop.table(frequencia_atividades) * 100

frequencia_nivh <- table(dados$NivH)
percentual_nivh <- prop.table(frequencia_nivh) * 100

cat("2) Frequência e percentual para Gender:\n")
print(frequencia_gender)
print(percentual_gender)
cat("\n")

cat("Frequência e percentual para Idade_G:\n")
print(frequencia_idade_g)
print(percentual_idade_g)
cat("\n")

cat("Frequência e percentual para Educ:\n")
print(frequencia_educ)
print(percentual_educ)
cat("\n")

cat("Frequência e percentual para Atividades:\n")
print(frequencia_atividades)
print(percentual_atividades)
cat("\n")

cat("Frequência e percentual para NivH:\n")
print(frequencia_nivh)
print(percentual_nivh)
cat("\n")

# Tests Pre - by group

# Shapiro-Wilk
dados %>% select(2:4) %>% group_by(NivH) %>% 
  doo(~mshapiro_test(.))


## Norm_univariada_by group
dados %>% group_by(NivH) %>% 
  shapiro_test(EGT, HPO)

## Abord1 Mahalanobis (outlier = p<0,001)

dados %>% dplyr::select(2:4) %>% group_by(NivH) %>% 
  doo(~mahalanobis_distance(.)) %>% 
  filter(is.outlier == TRUE)

## outliers UNIVARIADOS - por grupo:

boxplot(dados$HPO ~ dados$NivH)
boxplot(dados$EGT ~ dados$NivH)

dados %>% group_by(NivH) %>% 
  identify_outliers(HPO)

dados %>% group_by(NivH) %>% 
  identify_outliers(EGT)

# CORRELACAO Moderada ?

correlation <- cor(dados$EGT, dados$HPO)
cor_ci <- cor.test(dados$EGT, dados$HPO, conf.level = 0.95)

bootstrap_corr <- function(data, indices) {
  sampled_data <- data[indices, ]
  corr <- cor(sampled_data$EGT, sampled_data$HPO)
  return(corr)
}

set.seed(123)  
bootstrap_results <- boot(data = dados, statistic = bootstrap_corr, R = 1000)
bca_ci <- boot.ci(bootstrap_results, type = "bca")

print(paste("Correlação entre EGT e HPO:", round(correlation, 3)))
print(paste("Intervalo de Confiança (95%): [", round(cor_ci$conf.int[1], 3), ",", round(cor_ci$conf.int[2], 3), "]"))
print("Intervalo de Confiança BCa (95%):")
print(bca_ci)

#fim da correlação


# Homogeneidade
## Alfa: 0,001
box_m(dados[,c("HPO", "EGT")], dados$NivH)

## Levene
leveneTest(HPO ~ NivH, dados, center = mean)
leveneTest(EGT ~ NivH, dados, center = mean)

## Matriz de correlacao abordagem2
matriz <- cor(dados[,3:4])
View(matriz)

# Relacao linear var_depend_group
pairs(dados[,3:4], pch = 19,
      col = dados$NivH)

dplyr::select(dados, 2:4)

graf <- dados %>% dplyr::select(2:4) %>% group_by(NivH) %>% 
  doo(~ggpairs(.), result = "plots")

graf$plots

graf$plots[which=1]
graf$plots[which=2]
graf$plots[which=3]


# MANOVA Model

modelo <- manova(cbind(EGT, HPO) ~ NivH, data = dados)

## Results

options(scipen = 999)

summary(modelo, test = "Wilks")

## ANOVA univariada:

summary.aov(modelo)


# Estimated Marginal Means

dados %>% emmeans_test(HPO ~ NivH, p.adjust.method = "bonferroni")

dados %>% emmeans_test(EGT ~ NivH, p.adjust.method = "bonferroni")


# Analise descritiva

ggplot(dados, aes(x = EGT, y = HPO, group = NivH, color = NivH)) +
  geom_point()

dados %>% group_by(NivH) %>% 
  get_summary_stats(HPO, EGT, type = "mean_sd")

# CI MÉDIAS Hpo com Bca

n_iterations <- 1000

# Função CI BCa
calc_ci_bca <- function(data) {
  result <- boot(data$HPO, statistic = function(x, indices) mean(x[indices]), R = n_iterations)
  ci <- boot.ci(result, type = "bca")
  return(ci$bca)
}

results_bca <- by(dados, dados$NivH, calc_ci_bca)

# Resultados
results_bca

# ic MÉDIAS EGT com Bca ******
n_iterations <- 1000

calc_ci_bca <- function(data) {
  result <- boot(data$EGT, statistic = function(x, indices) mean(x[indices]), R = n_iterations)
  ci <- boot.ci(result, type = "bca")
  return(ci$bca)
}

results_bca <- by(dados, dados$NivH, calc_ci_bca)

# ******* RESULTS
results_bca